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ABSTRACT 


Fourier’s  approach  to  the  problem  of  heat  conduction 
served  as  the  introduction  to  Fourier  series  and  Fourier 
integral.  This  led  later  to  the  discrete  Fourier  transform 
and  recently  to  the  fast  Fourier  transform.  The  fast 
Fourier  transform,  serving  as  an  expeditious  means  to 
pass  from  functions  of  time  to  functions  of  frequency, 
contributes  to  the  analysis  of  time-invariant  systems  for 
deterministic  as  well  as  stochastic  situations. 


ADMINISTRATIVE  INFORMATION 

The  material  of  this  report  was  prepared  as  an  invited  tutorial 
talk  for  presentation  at  the  Eighty- Fifth  Meeting  of  the  Acoustical  Society 
of  America,  held  10-13  April  1973,  in  Boston,  Massachusetts.  Funding 
was  provided  under  Task  Area  SF  53  532,  Work  Unit  1-1820-002. 
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HISTORICAL  REMARKS  ABOUT  THE 
FOURIER  TRANSFORM 


The  development  of  the  Fourier  transform  goes  back  about  150  years 
to  Jean  Baptiste  Joseph  Fourier  who  was  engaged  at  the  time  in  trying 
to  solve  the  problem  of  heat  conduction.  The  partial  differential  equation 
of  heat  conduction  is  solved,  for  certain  regions,  by  separation  of 
variables.  This  leads  among  other  things,  to  an  ordinary  differential 
equation  which  calls  for  a  function  proportional  to  its  second  derivative, 
the  factor  of  proportionality  being  a  negative  number.  Functions 

I 

satisfying  such  differential  equations  are,  of  course,  sinusoidal.  To 
adjust  solutions  to  the  specific  conditions  of  a  particular  heat  problem, 
Fourier  developed  his  famous  theorem  that  any  quite  arbitrary  and 
irregular  function,  either  periodic  or  given  only  in  a  finite  interval,  can 
be  expanded  into  an  infinite  sum  of  cosines  and  sines  whose  arguments 
are  integer  multiples  of  one  basic  argument.  If  a  function  g(t)  is  periodic 
with  period  T,  or  is  given  only  for  the  interval  0  <  t  <  T,  then  g(t)  can 
be  written  as 


g(t)  = 


£ 


n=-co 


a 


n 


nt 


where 


a 

n 


1_ 

T 


dt 


This  relationship  is  expressed  here  by  means  of  exponentials,  making 
use  of  the  well-known  relation  e  =  cos  b  +  i  sin  b,  rather  than  by  means 
of  sines  and  cosines.  The  question  of  precisely  describing  those  functions 
which  equal  a  convergent  Fourier  series  has  occupied  mathematicians 
ever  since  the  time  of  Fourier.  Fourier’s  method  of  expanding  periodic 
functions  into  series  was  extended  to  functions  which  are  not  periodic. 

This  led  to  the  Fourier  integral,  customarily  written  as 
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g(t)  =  J  G(f)  e2?rift  dt 

-05 

where 

G(f)  =  Jg(t)  e"2?rift  dt  . 

-05 

00 

These  formulas  hold*  of  course,  only  if  the  improper  integral  J  |g(t)|  dt 
is  convergent.  They  associate  a  function  of  time,  t,  with  a 
function  of  frequency,  f.  Until  recently  the  term  Fourier  transform  was 
reserved  for  these  formulas. 

If  we  are  interested  in  numlerical  results  we  usually  deal  with  a 
situation  in  which  a  function  is  given  at  N, usually  equidistant, points: 
g(0),  g(l),  g(2),  „ .  „ ,  g(N-2),  g(N-l),  and  we  want  to  find  N quantities 
G(0),  G(l),  . ..,  G(N-l)  such  that 


g(k)  = 


N-l 

L 

q-0 


2iri 
G{q)  e  N 


kq 


Orthogonality  relations  permit  this  system  of  N  equations  in  N  unknowns 
to  be  readily  solved  for  the  G(q).  Thus 


G(q)  = 


1 

N 


N-l 

L  g(k)  e 
k=0 


2iri 

N 


kq 


The  computation  of  the  G(q),  which  we  call  the  discrete  Fourier  transform, 
from  given  g(k)  is  a  very  straightforward  process  involving  N 

multiplications  and  additions  to  find  one  G(q).  To  find  all  N  values  of 

2  2 
G(q),  N  operations  are  necessary;  since  N  for  large  N  is  a  very  large 

number,  a  shorter  method  was  needed  even  for  high-speed  computers, 
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and  the  fast  Fourier  transform  was  developed. 


THE  FAST  FOURIER  TRANSFORM 

A  substantial  saving  of  time  can  be  achieved  if  the  integer  N  has 
many  factors.  There  are  various  ways  of  describing  the  time  saving 

process  of  the  fast  Fourier  transform  (FFT),  one  of  which  follows: 

6 

Let  N  be  a  power  of  2,  for  example ,  N  =  64  =  2  .  Then  the  original 
problem  is  equivalent  to  multiplying  a  64  dimensional  vector  by  a  64  by  64 
matrix.  None  of  the  elements  of  the  matrix  are  zero,  since  they  are  all 
of  the  form  e^  which  has  an  absolute  value  of  1.  It  can,  however,  be 
shown  that  this  matrix  can  be  factored  into  two  matrices,  one  having 
32  non-zeros  per  row  and  per  column,  and  the  other  having  two  non-zeros 
per  row  and  per  column.  We  can  then  further  factor  the  first  matrix, 
the  one  with  the  32  non-zeros  in  each  row  and  each  column,  and  continue 
the  process  until  we  obtain  six  matrices  which  have  only  two  non-zero 
elements  in  each  row  and  in  each  column.  Thus,  instead  of  needing 
64  =  4096  operations,  we  now  need  only  2  x  64  x  6  or  768  operations, 
a  substantial  reduction  in  the  number  of  operations  required.  Note  that 
reference  to  the  matrices  was  made  merely  to  illustrate  how  the  FFT 
saves  time.  There  is  no  intention  to  suggest  that  the  FFT  should  be 
carried  out  by  actually  multiplying  the  matrices  described. 


Cooley,  J.W.  andJ.W.  Tuckey,  "An  algorithm  for  the  machine 
calculation  of  complex  Fourier  Series,  "  Math.  Comp.,  vol.  10,  pp.  297- 
301,  April  1965. 
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APPLICATIONS 


The  following  applications  have  been  chosen  to  illustrate  the  use  of 
the  Fourier  transform  in  general  and  the  fast  Fourier  transform  in 
particular.  No  effort  has  been  made  to  select  the  most  important  or  the 
most  recent  applications.  Note  that  the  fast  Fourier  transform  may  be 
used  to  advantage  in  any  large  problem  in  which  the  Fourier  transform  is 
being  used. 

LINEAR  DIFFERENTIAL  EQUATIONS  WITH  CONSTANT  COEFFICIENTS 

One  of  the  most  important  applications  of  the  Fourier  transform, 
fast  or  otherwise,  has  to  do  with  the  solution  of  linear  ordinary  differential 
equations  with  constant  coefficients.  For  instance  we  study 


S  b  g(m)(t)  =  p(t) 

0  111 


where  p(t)  is  a  given  function  and  bm  are  given  coefficients,  and  where 


g(m)(t) . 


dt 


The  differential  equation  is  to  be  solved  for  the 


unknown  function  g(t) .  Applying  the  Fourier  transform  to  both  sides  of 
the  differential  equation  changes  the  left-hand  side  into  a  product  of  the 
transform  G(f)  and  a  function  of  the  frequency  f,  such  that  we  may  write 
G(f)  =  H(f)  •  P(f)  where  P(f)  is  the  transform  of  p(t).  From  G(f)  we  then 
go  to  g(t)  by  taking  the  inverse  Fourier  transform.  To  use  this  method, 
it  is,  of  course,  important  to  be  able  to  move  easily  from  function 
of  time  to  function  of  frequency,  and  vice  versa,  and  it  is  this  facility 
that  the  FFT  provides. 
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LINEAR  TIME-INVARIANT  SYSTEMS 


Often  we  are  in  a  worse  predicament  in  which  we  are  to  solve  the 
equation 

£  bm  =  P(0 

without  actually  being  given  the  coefficients  b  .  We  are  given  only  a 
physical  situation  which  takes  input  functions  and  furnishes  output 
functions.  Such  a  physical  system  which  is  described  by  a  linear  differen¬ 
tial  equation  with  constant  coefficients  is  called  a  linear  time -invariant 
system.  The  problem  is  to  find  g(t)  if  p(t)  is  given.  We  have  to  find 
the  function  H(f)  without  using  the  coefficients  b  .  By  assuming  that 
we  have  before  us  a  physical  system,  we  make  an  experiment  with  an 
input  which  is  sinusoidal  with  a  particular  frequency,  say  f^.  Then  a 
simple  analysis  shows  that  H^q)  becomes  the  ratio  of  output  to  input. 

The  disadvantage  of  this  method  is  that  it  gives  the  function  H(f)  for 
only  one  particular  f;  so  then  for  each  f  desired,  another  experiment 
must  be  made. 

A  more  expeditious  way  of  finding  H(f)  would  be  to  choose  p(t)  such 
that  P(f)  =  1,  so  that  the  corresponding  G(f)  will  be  the  function  H(f) 
which  governs  the  system.  But  it  is  shown  in  the  standard  books  on 
Fourier  transforms  that  the  inverse  Fourier  transform  of  the  constant  1 
is  the  Dirac  function.  If  we  find  the  response  of  our  system  to  a  Dirac 
function,  we  will  have  a  function  h(t),  which  we  call  the  impulse  response 
whose  Fourier  transform  is  H(f),  the  so-called  transfer  function.  Now 
if  we  want  to  find  the  g(t)  that  goes  with  a  given  p(t),  we  first  find 
P(f),  the  Fourier  transform  of  p(t).  Then  G(f)  =  H(f)  •  P(f),  and  g(t) 
is  the  inverse  Fourier  transform  of  G(f). 
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A  practical  application  of  this  method,  described  to  me  by  a  colleague, 

2 

Theo  Kooij,  concerns  the  prediction  of  the  echo  from  an  underwater 

object  for  various  sonar  signals.-' .  Assuming  that  the  echo  can  be 

considered  as  output  of  a  linear  time -invariant  system,  we  find  the 

response  to  a  Dirac  function  by  finding  the  echo  to  the  explosion  of  an 

explosive  charge.  This  echo  is  the  function  h(t)  associated  with  the 

system  (or  more  exactly  a  constant  multiple  thereof).  We  digitize  the 

hydrophone  output  and  find  the  fast  Fourier  transform,  which  gives  us 

H(f).  For  any  sonar  signal  whose  echo  we  want  to  find, -weifirSt  find  the 

Fourier  transform,  which  can  be  done  explicitly  in  many  cases,  and  then 

multiply  it  by  H(f)  and  take  the  fast  inverse  Fourier  transform  of  the  product. 

Thus  the  echo  for  various  sonar  signals  can  readily  be  found. 

Another  example  of  a  time -invariant  linear  system  is  in  connection 

with  the  response  of  a  non-uniform  beam  to  a  transient  forcing  function. 

Such  a  Situation  is  described  by  Francis  Henderson  in  his  report, 

"Transient  Response  Calculation  in  the  Frequency  Domain  with  General 

g 

Bending  Response  Program.  "  In  this  case  the  unknown  function  is  either 
the  displacement  or  the  moment,  each  a  function  of  two  variables,  time  t 
and  longitudinal  coordinate  x.  Here  we  have  a  differential  equation 
which,  after  applying  the  Fourier  transform  with  respect  to  the  variable  t, 
still  contains  derivatives  with  respect  to  the  variable  x.  A  difference 
equation  scheme  is  used  by  dividing  the  x  interval  into,  say,  M  sub¬ 
intervals.  The  differential  equation  then  leads  to  a  system  of  M  linear 
equations  which  contain  the  frequency  f  as  a  parameter.  This  system 
must  be  solved  for  a  large  number  of  values  of  the  frequency  f.  These 


Personal  Communication 
3 

Henderson,  F. ,  "Transient  Response  Calculation  in  the  Frequency 
Domain  with  General  Bending  Response  Program,  "  Naval  Ship  Research 
and  Development  Center  Report  1613,  February  1971. 
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solutions  are  obtained  by  utilizing  the  fact  that  the  matrix  in  question  has 
zeros  in  most  positions  except  in  the  vicinity  of  the  main  diagonal.  After 
having  obtained  the  solution  as  functions  of  frequency,  we  then  obtain  the 
solution  as  a  function  of  time  by  means  of  the  inverse  fast  Fourier 
transform. 

STOCHASTIC  PROCESSES 

In  the  foregoing  cases,  a  linear  time -invariant  device  connected  an 

input  and  an  output  function.  There  are,  however,  many  physical  situations 

in  which  it  is  more  meaningful  to  consider  the  input  and  the  output,  not 

as  deterministic  functions  of  time,  but  as  stochastic  processes.  For 

instance,  we  might  be  interested  in  studying  the  interrelation  between  the 

heights  of  the  waves  in  the  ocean  and  the  response  of  a  ship.  One  such 

4 

study  has  been  made  by  M.  K.  Ochi  and  W.  E.  Bolton,  in  which  they 
investigate  the  interrelation  between  the  input  stochastic  process,  given 
by  the  wave  height  at  a  chosen  point  in  the  ocean,  and  the  output  stochastic 
process,  the  response  of  the  ship.  The  connection  between  input  and 
output  stochastic  processes  can  again  be  described  by  means  of  the 
impulse  response  h(t)  or  the  transfer  function  H(f),  its  Fourier  transform, 
but  the  formulas  look  a  bit  different.  Here,  S(f),  the  so-called  power 
spectral  density  —  often  called  simply  the  power  spectrum  —  of  the  input 
stochastic  process  is  related  to  U(f),  the  power  spectrum  of  the  output 
stochastic  process,  by  means  of 

U(f)  =  |  H(f)  | 2  S(f) . 

The  power  spectrum  S(f),  and  likewise  U(f),  is  a  function  of  frequency 
characterized  by  the  property  that  the  integral  over  a  frequency  interval 


Ochi,  M.  K.  andW.E.  Bolton,  "Statistics  for  Prediction  of  Ship 
Performance  in  a  Seaway,  "  International  Shipbuilding  Progress,  Vol.  20, 
No.  222  (1973). 
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(fpf2)  namely  J  S(f)  df  represents  the  contribution  to  the  variance 
fl 

from  that  frequency  interval.  (This  definition  by  an  integral  explains 

why  the  longer  name  "power  spectral  density”  is  somewhat  more  precise 

than  "power  spectrum". )  The  power  spectrum  can  be  found  as  the 

Fourier  transform  of  the  autocorrelation  which  in  turn  is  found  from  an 

observed  particular  realization  of  the  stochastic  process.  An  alternate 

way  of  computing  the  power  spectrum,  and  here  we  refere  to  the  discrete 

case  for  N  given  points,  is  to  first  find  the  Fourier  transform  of  the 

given  data,  and  to  then  find  the  squares  of  the  absolute  values.  Before 

the  introduction  of  the  FFT  this  latter  method  was  found  to  be  much 

more  time  consuming  than  finding  the  autocorrelation  for  a  certain  number 

of  lags  and  then  its  Fourier  transform.  However,  the  FFT  now  makes  this 

alternate  method,  which  does  not  use  the  autocorrelation,  much  faster 

than  the  method  which  uses  autocorrelation. 

Returning  to  our  discussion  of  the  interrelation  of  two  stochastic 

processes  (the  wave  height  in  the  ocean  and  the  response  of  a  ship),  let  us 

consider  a  possible  way  of  finding  the  function  H(f),  or  the  square 
2 

|  H(f)  |  ,  known  among  naval  engineers  as  the  response  amplitude  operator. 
If  we  create  waves  of  a  certain  frequency  around  a  ship  model,  we  can 
observe  the  response  of  the  model,  which  may  be  a  displacement,  a 
force,  a  moment,  whatever  quantity  we  want  to  study.  Then  H(f),  the 
transfer  function  at  a  particular  frequency,  equals  the  ratio  of  output 
to  input  for  a  sinusoidal  wave.  Once  the  function  H(f)  is  determined,  or 
more  precisely,  once  the  functions  H(f)  with  respect  to  each  interesting 
response  are  determined,  the  corresponding  power  spectra  of  the 
responses  are  readily  found  for  any  given  power  spectrum  of  the  waves 
of  the  ocean. 
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SIGNAL  ANALYSIS 


We  have  thus  far  given  examples  of  finding  the  relation  between  the 
input  time  function  and  the  output  time  function  by  first  finding  the  input 
as  a  function  of  frequency  and  then  by  a  mere  multiplication  finding  the 
output  first  as  function  of  frequency  and  then,  by  inverse  Fourier 
transform,  as  function  of  time.  Sometimes,  however,  we  have  problems 
which  do  not  require  that  all  of  these  steps  be  taken.  For  instance, 
suppose  we  are  listening  to  underwater  sound  in  order  to  find  out  whether 
there  are  any  interesting  sound  sources  in  the  water  which  we  would  like 
to  identify.  We  can  consider  the  underwater  sound  as  a  stochastic  process 
and  use  the  time-series  obtained  by  digitizing  the  output  of  a  hydrophone 
to  compute  the  power  spectrum  of  the  underwater  sound.  We  need  not 
subject  the  input  stochastic  process  to  a  linear  time-independent  operation, 
but  can  be  satisfied  with  studying  the  power  spectrum  of  the  underwater 
sound.  We  observe  the  characteristic  features  of  the  power  spectrum,  the 
location  of  its  peaks  and  the  like,  and  use  this  observation  to  compare  it 
with  previously  obtained  power  spectra  and  store  its  essential  characteristics 
for  later  comparisons.  In  this  problem,  we  need  only  to  perform  a  direct 
Fourier  transform,  the  inverse  Fourier  transform  is  not  needed. 

A  quite  similar  application  has  to  do  with  listening  to  one  of  our  own 
vessels.  We  again  find  the  power  spectrum  and  identify  from  it 
objectionable  sources  of  sound  on  the  vessel  with  the  view  of  making  the 
vessel  more  nearly  silent.  In  this  application,  many  of  the  computer 
programs  developed  for  detection  of  objects  by  underwater  sound  are 
used. 

TURBULENCE  STUDIES 

One  of  the  early  applications  of  the  power  spectrum  was  in  the  study 
of  turbulence.  We  are  concerned  with  a  fluid'  motion  for  which  it  is 
meaningful  to  consider  the  velocity  as  a  stochastic  process.  To  get  better 
insight  into  the  phenomenon  we  use  the  FFT  to  determine  the  power 
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spectrum.  In  calling  the  fluid  motion  a  stochastic  process  we  need  to  give 
some  further  explanation.  Thus  far  we  have  referred  to  stochastic 
processes  which  depend  on  time  and  on  chance.  A  phenomenon  depending 
on  chance  and  on  a  coordinate  x  may  just  as  well  be  called  a  stochastic 
process.  In  general,  we  may  study  stochastic  processes  that  depend 
on  chance,  on  time,  and  on  the  three  coordinates  x,  y,  z.  In  the  case 
of  a  function  of  time  we  are  led  to  a  Fourier  transform  depending  on 
frequency,  where  frequency  has  the  dimension  of  time  to  the  power  minus 
one;  frequency  gives  the  number  of  waves  per  unit  of  time.  A  function  of 
x,  of  length,  yields  a  Fourier  transform  depending  on  inverse  length, 
which  may  be  called  spatial  frequency  or  wave  number,  giving  the  number 
of  waves  per  unit  length.  In  the  study  of  turbulence  we  deal  with  stochastic 
processes  depending  on  time  as  well  as  with  those  depending  on  space. 

PICTURE  PROCESSING 

Another  application  of  Fourier  analysis  to  functions  of  spatial 
coordinates  is  met  in  the  field  of  picture  processing.  If,  for  instance,  we 
want  to  remove  what  may  be  called  optical  noise  from  a  picture,  we 
can  often  do  so  by  filtering,  which  in  turn  is  done  by  first  finding  the 
Fourier  transform,  then  removing  components  of  undesirable  spatial 
frequencies  and  finally  taking  the  inverse  Fourier  transform. 

CONCLUDING  REMARKS 

In  conclusion,  attention  should  be  called  to  the  fact  that  the  applications 
of  the  FFT  are,  of  course,  also  applications  of  high-speed  computers.  It 
is  no  accident  that  the  FFT  was  developed  at  a  time  when  means  were 
available  to  attack  problems  involving  very  many  points,  in  situations 
where  the  FFT  could  show  its  vast  superiority  over  other  methods. 
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